---
title: "Fairwheather Cosmo: Main Analysis"
output:
  html_document:
    df_print: paged
  html_notebook: default
  pdf_document: default
---



```{r setup, warning=F, message=F}

# Load pacman
library(pacman)

# Use pacman to load libraries
p_load(tidyverse, haven, lmerTest, xtable, janitor, knitr, huxtable, marginaleffects,
       ggplot2, ggpubr, ggsci, countrycode, readxl, texreg, dotwhisker, ggrepel,
       here, psych, ggcorrplot, Hmisc, xtable)

# Read the WVS data
wvs_latam <- read_csv("wvs_latam.csv")

# Make the var names lowercase
wvs_latam <- clean_names(wvs_latam)


```


## Alpha of Immigration Index 
 
```{r alpha}

alpha_imm <- wvs_latam %>% 
  select(starts_with("imm"), -c(imm_index)) %>% 
  alpha(check.keys = T)

pluck(alpha_imm, "total") %>% 
  data.frame() %>% 
  select(raw_alpha)

```

## Plot Immigration Index by Country (Figure 1)


```{r plot-immig-index}


# Means by country 
mean_imm <- wvs_latam %>% 
  group_by(country_str) %>% 
  summarise(mean_index = mean(imm_index, na.rm =T))

# Join means to the main data
wvs_latam <- left_join(wvs_latam, mean_imm)

wvs_latam <- wvs_latam %>% 
  drop_na(country_str)

wvs_latam$country_str <- fct_reorder(wvs_latam$country_str, wvs_latam$mean_index)

ggplot(wvs_latam, aes(x = imm_index, fill = "density")) + 
  geom_density(alpha = .25) +
  facet_wrap(~country_str) +
  geom_vline(aes(xintercept = mean_index), 
             linetype = "dotted", color = "red")  + 
  geom_text(aes(x = mean_index, y = 0.03, label = round(mean_index, 1)), 
            color = "red", hjust = 0, vjust = 3) +
  labs(x = "Immigration Index", y = "Density") +
  scale_y_continuous(limits = c(0, 0.03)) +
  scale_x_continuous(breaks = round(seq(20, 100, by =20),1), limits = c(10, 100)) +
  scale_fill_manual(values = "lightblue", guide = "none") + 
  theme_pubr(base_size = 14) +
  theme(axis.title = element_text(face = "bold"))

ggsave(file = here("figures/immig_by_country.pdf"),
       height = 6,
       width = 8)

  
```


Immigrant v refugee measures:

```{r immig_v_migrant}

wvs_latam %>% 
  select(iso3c, refugee_change, mig_change) %>% 
  unique() 

wvs_latam %>% 
  select(refugee_change, mig_change) %>% 
  unique() %>% 
  cor()

```



# Effects by Country (Figures 2 and 3)

Let's run the model by country and plot the estimated effect of cosmopolitanism and national identity against the value of refugee change.

## Cosmopolitanism 

Manuscript Figure 2:

```{r effects-cosmo, warning=FALSE}

# Prep the data 
cosmo_dat <- wvs_latam %>% 
  select(iso3c, imm_index, close_world, refugee_change) %>%  # select vars
  na.omit() %>% # drop missings 
  group_by(iso3c, refugee_change) %>%  # group by country
  nest() # nest the data

# Function for our model
cosmo_mod <- function(df) {
  lm(imm_index ~ close_world, data = df) }

# Function for CIs 
ci <- function(mod) {
  confint(object = mod, parm = "close_world")
}


# Run the model on each country
cosmo_dat <- cosmo_dat %>% 
  mutate(model = map(data, cosmo_mod),  # run the model by country
         coeffs = map(model, coefficients),# extract the coefficents 
         conf = map(model, ci))  # get confidence intervals

# Grab the point estimate
cosmo_dat$cosmo <-  sapply(cosmo_dat$coeffs, "[[", 2)

# Grab the lower CI
cosmo_dat$lower <- sapply(cosmo_dat$conf, function(matrix) {
  matrix[, 1] })

# Grab the upper CI
cosmo_dat$upper <- sapply(cosmo_dat$conf, function(matrix) {
  matrix[, 2] })


#Plot effect of cosmopolitanism against refugee change
ggplot(cosmo_dat, aes(x = refugee_change,y = cosmo, label = iso3c)) + 
  geom_point() +
  geom_smooth(method = "lm", se = F, linetype = "dotted") +
  geom_linerange(aes(ymax = upper, ymin = lower), alpha = .30, size = .60) +
  geom_text_repel(min.segment.length = 0) +
  labs(x = "% Change in Refugee Population", 
       y = expression(bold(paste(hat(beta)," on Cosmopolitanism")))) +
  theme_pubr(base_size = 14) +
  theme(axis.title = element_text(face = "bold"))

ggsave(file = here("figures/cosmo_by_country.pdf"),
       height = 6,
       width = 8)


```


## National Pride 

Manuscript Figure 3:

```{r effects-nationalism, warning=FALSE}

# Prep the data 
nat_dat <- wvs_latam %>% 
  select(iso3c, imm_index, natpride, refugee_change) %>%  # select vars
  na.omit() %>% # drop missings 
  group_by(iso3c, refugee_change) %>%  # group by country
  nest() # nest the data

# Function for our model
nat_mod <- function(df) {
  lm(imm_index ~ natpride, data = df) }

# Function for CIs 
ci_nat <- function(mod) {
  confint(object = mod, parm = "natpride")
}


# Run the model on each country
nat_dat <- nat_dat %>% 
  mutate(model = map(data, nat_mod),  # run the model by country
         coeffs = map(model, coefficients), # extract the coefficents 
         conf = map(model, ci_nat)) #


# Run the model on each country
nat_dat <- nat_dat %>% 
  mutate(model = map(data, nat_mod),  # run the model by country
         coeffs = map(model, coefficients)) # extract the coefficents 


# Grab the coefficient on national pride 
nat_dat$nat <-  sapply(nat_dat$coeffs, "[[", 2)

# Grab the lower CI
nat_dat$lower <- sapply(nat_dat$conf, function(matrix) {
  matrix[, 1] })

# Grab the upper CI
nat_dat$upper <- sapply(nat_dat$conf, function(matrix) {
  matrix[, 2] })

ggplot(nat_dat, aes(x = refugee_change,y = nat, label = iso3c)) + 
  geom_point() +
  geom_smooth(method = "lm", se = F, linetype = "dotted") +
  geom_linerange(aes(ymax = upper, ymin = lower), alpha = .35, size = .65) +
  geom_text_repel(min.segment.length = 0) +
  labs(x = "% Change in Refugee Population", 
       y = expression(bold(paste(hat(beta)," on National Identity")))) +
  theme_pubr(base_size = 14) +
  theme(axis.title = element_text(face = "bold"))

ggsave(file = here("figures/nat_by_country.pdf"),
       height = 6,
       width = 8)


```


##  Descriptives for SI

```{r descriptives}

# Quick function to drop na's and take mean
na_mean <- function(x) {
  return(mean(na.omit(x)))
}


# Number of obs and mean imm_index
wvs_desc1 <- wvs_latam %>% 
  group_by(iso3c) %>%  # by country
  summarise(
    obs = n(), # obs 
    mean_imm = na_mean(imm_index), # mean immigration index
    mean_cosmo = na_mean(close_world), # mean cosmo
    mean_nat = na_mean(natpride), # mean nationalism
    ) %>%  
  select(iso3c, N = obs, "Imm. Index" = mean_imm,
         "x-bar cosmo" = mean_cosmo, "x-bar nat" = mean_nat) 

# Year of the wvs survey
wvs_desc2 <- wvs_latam %>% 
  select(iso3c, year, "% Change Refugees" = refugee_change) %>%  # select country and year
  unique.data.frame() # grab unique country year combos

# Combine the descriptives
wvs_desc3 <- full_join(wvs_desc1, wvs_desc2, by = "iso3c")  %>% 
  select(iso3c, year, everything())

# Add the estimated betas on cosmo
wvs_desc4 <- full_join(wvs_desc3,  cosmo_dat[c("iso3c", "cosmo")])

# Add the estimated betas on nat
wvs_desc5 <- full_join(wvs_desc4, nat_dat[c("iso3c", "nat")])

# Table 
kable(wvs_desc5, digits = 2)

#Out to latex
print(xtable(wvs_desc5), include.rownames = FALSE)

```


# Main WVS Models 

## Recoding/Merging for Models 

```{r shaping}

# Make Education a Factor
wvs_latam <- wvs_latam %>% 
  mutate(edu_cat = case_match(edu,
                              1:2 ~ "Primary", 
                              1:4 ~ "Secondary",
                              5:8 ~ "Post-Secondary"))

wvs_latam$edu_cat <- factor(wvs_latam$edu_cat, levels = c("Primary", 
                                              "Secondary", 
                                              "Post-Secondary"))


# Log the Income Level Variable 
wvs_latam$income_log <- log(wvs_latam$income_lvl)

# Log Age 
wvs_latam$ln_age <- log(wvs_latam$age)

# Make a couple of variables factors
wvs_latam$unemp <- factor(wvs_latam$unemp) #unemployment
wvs_latam$female <- factor(wvs_latam$female) #female
wvs_latam$urban <-  factor(wvs_latam$urban) # urban rural


```


## Models 

- Multilevel models with country-level random effects. 
- First estimate main effect of cosmopolitanism and national pride, then estimate interaction with change in number of refugees. 

### Manuscript Table 1 

```{r estimate-mods}

# Main Effect of Cosmo
m1 <- lmer(imm_index ~ close_world + edu_cat + trust_generalized + income_log + urban +
             ln_age + female + unemp + refugee_change + ln_gdppc + (1| iso3c), data = wvs_latam) 


# Predictions: min to max
cosmo_m1 <-  predictions(m1, newdata = datagrid(close_world = c(0, 3))) %>% 
  select(estimate, close_world)

cosmo_m1[2,1] - cosmo_m1[1,1] 


# Conditional Effect of Cosmopolitanism 
m2 <- lmer(imm_index ~ close_world * refugee_change + edu_cat + trust_generalized + 
             income_log + urban + unemp + ln_age +  female + unemp + ln_gdppc +  (1| iso3c), data = wvs_latam) 

summary(m2)

# Main Effect of National Pride
m3 <- lmer(imm_index ~ natpride + edu_cat + trust_generalized + income_log + ln_age + female + 
             unemp + urban + refugee_change + ln_gdppc + (1| iso3c), data = wvs_latam) 

  

nat_m3 <- predictions(m3, newdata = datagrid(natpride = c(0, 4))) %>% 
  select(estimate, natpride) 
  

nat_m3[1,1] - nat_m3[2,1]


# Conditional Effect of National Pride
m4 <- lmer(imm_index ~ natpride * refugee_change +  edu_cat + trust_generalized + income_log
           + urban + ln_age + female + unemp + ln_gdppc + age + female +  (1| iso3c), data = wvs_latam)

summary(m4)

# List of our models 
mods <- list(m1, m2, m3, m4)

# Coefficient Labels 
coef_labs <- list(
                  "close_world" = "Cosmopolitanism", 
                   "close_world:refugee_change" = "Cosmopolitanism $\\times$ Refugee Change",
                  "natpride" = "National Pride",
                  "natpride:refugee_change" = "National Pride $\\times$ Refugee Change",
                  "trust_generalized" = "Interpersonal Trust",
                  "edu_catSecondary" = "Secondary Edu", 
                  "edu_catPost-Secondary" = "Post-Secondary Edu", 
                  "ln_age" = "Age (logged)",
                  "urban1" = "Urban",
                  "unemp1" = "Unemployed",
                  "female1" = "Female",
                  "income_log" = "Individual Income", 
                  "refugee_change" = "Refugee Change", 
                   "ln_gdppc" = "Log GDP/Capita")

# Print Models
screenreg(mods,
          stars = c(0.1, 0.05, 0.01),
          custom.coef.map = coef_labs)


# Caption 
cap <- "Multilevel models of pro-immigrant attitudes. Random intercepts at the country level. Higher values of the DV indicate more pro-immigrant attitudes."

#Out to Latex
texreg(list(m1, m2, m3, m4),
       custom.coef.map = coef_labs,
       stars = c(0.1, 0.05, 0.01),
       caption = cap,
       file = "tables/mods_z.tex",
       label = "tab_main_table")


```


##  Marginal Effects (Figures 4 and 5)

```{r plot-margins}

# Marginal Effect of Cosmopolitanism 
me_cosmo <- plot_slopes(model= m2, 
            variables = "close_world", 
            condition = "refugee_change", 
            rug = T) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "% Change in Refugee Population", y = "Effect of Cosmopolitanism") +
  theme_pubr(base_size = 14) +
  theme(axis.title = element_text(face = "bold"))


me_cosmo

ggsave(file = here("figures/me_cosmo.pdf"),
       height = 6,
       width = 8)

# Grab the marginal effects for discussion
me_cosmo_dat <- me_cosmo[["plot_env"]][["dat"]]



# Marginal Effect of National Pride
me_nat <- plot_slopes(model= m4, 
            variables = "natpride", 
            condition = "refugee_change", 
            rug = T) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "% Change in Refugee Population", y = "Effect of Nationalism") +
  theme_pubr(base_size = 14) +
  theme(axis.title = element_text(face = "bold"))

me_nat

ggsave(file = here("figures/me_nat.pdf"),
       height = 6,
       width = 8)


me_nat_dat <- me_nat[["plot_env"]][["dat"]]


```


